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We report on a calculation of the B*Bk coupling in three-flavour lattice QCD. This coupling, defined from 
the strong-interaction matrix element {Bk\B*), is related to the leading order low-energy constant in heavy 
meson chiral perturbation theory (HM;t;PT). We carry out our calculation directly at the fc-quark mass using 
a non-perturbatively tuned clover action that controls discretization effects of order \pa\ and (ma)” for all n. 

Our analysis is performed on RBC/UKQCD gauge configurations using domain-wall fermions and the Iwasaki 
gauge action at two lattice spacings of a^* = 1.729(25)GeV, a^* = 2.281(28)GeV, and unitary pion masses 
down to 290MeV. We achieve good statistical precision and control all systematic uncertainties, giving a final 
result for the coupling = 0.56(3)stat(7)sys in the continuum and at the physical light-quark masses. This is the 
first calculation performed directly at the physical ii-quark mass and lies in the region one would expect from 
carrying out an interpolation between previous results at the charm mass and at the static point. 


I. INTRODUCTION 

The power of lattice QCD in probing the Standard Model 
and uncovering evidence for new physics lies predominantly 
in the flavour sector. To constrain the Cabibbo-Kobayashi- 
Maskawa (CKM) unitarity triangle [1-3] requires many in¬ 
puts that must be evaluated non-perturbatively, particularly in 
the B-meson sector. For instance, an important constraint on 
the apex of the CKM unitarity triangle comes from neutral B- 
meson mixing, which gives information on the ratio of CKM 
elements IVfiP/IVrdP- Accessing these CKM elements from 
the experimental data requires knowledge of the B-meson de¬ 
cay constant and bag parameter, or alternatively the SU(3) 
breaking ratio 

( 1 ) 

Lattice calculations of the decay constants fpj and are 
also necessary inputs for the Standard Model predictions of 
BR(B —?> tv) and BR(Bi —respectively, while lattice 
calculations of the B -A nlv form factor allow a determination 
of the CKM matrix element \Vuh\. For both semileptonic form 
factors and mixing matrix elements, the precision of lattice 
calculations lags behind experiment. The experimental mea¬ 
surements will continue to improve with the large data sets 
available at Belle II and from LHCb. Therefore it is essen¬ 
tial to reduce further the theoretical uncertainties in the non- 
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perturbative hadronic parameters in order to maximise the sci¬ 
entific impact of current and future B-physics experiments. 

A major source of uncertainty in all previous lattice calcula¬ 
tions is from practical difficulties simulating at physical light- 
quark masses. Theoretical insight from HM;kPT can guide ex¬ 
trapolations down to the physical point, but lack of knowledge 
of the low-energy constants (LECs) of the theory introduces 
uncertainties. For example, at next-to-leading order (NLO) in 
HMjPT and lowest order in the heavy-quark expansion the 
logarithmic dependence of and Bp^ on the light-quark (or 
equivalently, pion) mass is given by [4, 5] 

BB,^B[\-]^^^^Ml\og{Mllf^)^+---, (3) 

where g is the leading-order LEC. The strong-interaction ma¬ 
trix element (B7r|B*) is used to determine a coupling gt,, 
which would become g in the static limit of an infinitely heavy 
b quark. At the order used above in Eqs. (2) and (3), we are 
free to use gh in place of g; differences between the two are of 
order 1 /mi,. 

In this paper we perform the first calculation of the coupling 
gh directly at the ( 7 -quark mass. Previous determinations of the 
coupling have been hindered by the difficulties of simulating 
heavy quarks on the lattice. Lattice calculations have been 
performed for g^, the analogous coupling for D-mesons [6- 
9], and for g itself [7, 10-13]. Having a reliable theoretical 
calculation of the coupling for the B system is important since 
this coupling cannot be accessed directly through experiment. 
The strong coupling gD*D 7 i has been measured by the CLEO 
collaboration [14, 15] and more recently by BaBar [16, 17], 
but with B-mesons there is not enough phase-space for the 
B* ~A Bn decay to occur. Model estimates exist for gy, includ¬ 
ing from QCD sum rules [18-21] and non-relativistic quark 
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models [22]. 

The rest of this paper is organised as follows. In section II 
we briefly review the framework of HMjPT, show how gh 
enters and present the method for extracting gh from lattice 
matrix-element calculations. Section III describes the param¬ 
eters of the light-quark, gluon, and heavy-quark actions used 
in the numerical calculation and presents the ratios of two- 
and three-point correlators used to obtain gh- In section IV 
we fit the correlator ratios to extract gh and then extrapolate 
these results to the continuum and physical quark masses us¬ 
ing SU(2) HMjPT. We estimate the systematic errors in gh in 
section V, discussing each source of uncertainty in a different 
subsection. We conclude in section VI by presenting our fi¬ 
nal results and error budget, and comparing our result to other 
similar calculations. 


II. HEAVY MESON CHIRAL PERTURBATION THEORY 


In the infinite heavy-quark mass limit the properties of 
heavy-light mesons become independent of the heavy quark’s 
spin and flavour quantum numbers. Combining this with the 
chiral symmetry present in the massless light-quark (m^ —0) 
limit of QCD provides the basis for heavy meson chiral per¬ 
turbation theory. This effective theory of QCD is a joint ex¬ 
pansion in powers of the inverse heavy-quark mass 1 /mg and 
the light-quark-mass iriq. 

In HMjPT the heavy-light pseudoscalar and vector 
mesons, P and P*, are combined in a covariant 4x4 matrix 
representation 

H=^^{P;7^^-P75)■ (4) 


If one includes three light dynamical quark flavors {u,d,s), 
this corresponds to SU(3) HM;|^PT with the usual octet of 
pseudo-Nambu-Goldstone bosons for the light pseudoscalars; 






C6' 


K 


-1— ttO J-1_ ri 


K+ 


(5) 




Because the strange-quark mass is almost thirty times larger 
than the average up-down quark mass, however, one can also 
treat the strange quark as heavy and include only the up- 
and down-quark dynamical degrees of freedom; this leads 
to SU(2) HMjPT (with the corresponding modification of 
^). At lowest order the interactions between the heavy and 
light mesons are determined by a Lagrangian with a single 
LEC [23, 24] 

■^HMxPT ~ S'^^i^aHh■^^^^7^l75)^ ( 6 ) 

where 

(7) 

and ^ = exp{i^//%)■ The roman indices run over light-quark 
flavour and the trace is over Dirac indices. We use a conven¬ 
tion where « 130MeV. 


The matrix element for the strong transition B* Bn is 
parametrised by gB*BK, 

{B{p)n{q)\B*{p,X)) =gB*BKq-£‘'^\p), (8) 

where q = p — p' and e^‘^\p) is the polarization vector for 
polarization state labelled by A. Evaluating the same matrix 
element at leading order in HM;kPT, 

{P{p')Tl{q)\P*{p,X)) = ^^gq-e^^\p), (9) 

JTC 

enables the determination of gh from 


gB*BK = (10) 

Jit 

with gh equal to g up to 1/m]] corrections. 

Performing a Lehmann-Symanzik-Zimmermann reduction 
and using the partially-conserved axial current relation for a 
soft pion, Eq. (8) becomes 

gB-B7tq-£^^Hp) = 

j [/[2 _„2 f 

iqg^^ j d^xe'‘>-^{B{p')\Anx)\B*{p,X)), ( 11 ) 

where = vri 7 ^ 75 Vf 2 is the light-quark axial-vector current. 
Using a form-factor decomposition of the matrix element 


{B{p'W\B*{p,X)) = 2MB*Ao{q^ 




2 ‘l 


+ {MB*+MB)Ai{q^) 

g(A) 




e(^). 


q 




2 ^ 2(71 


■q 


Mb* +Mb 


■7 


we see that &tq^ = 0 


gB^BK = 


2M b*Ao{0) 
fit 


( 12 ) 


(13) 


On the lattice, we cannot simulate exactly at q^ = 0 with¬ 
out using twisted boundary conditions. Eurthermore and from 
Eq. (11), we see that the form factor Aq contains a pole at the 
pion mass, so it will be difficult to do a controlled extrapola¬ 
tion to q^ = 0. However, the decomposition in Eq. (12) must 
be free of unphysical poles, which allows us to obtain the re¬ 
lation 


gB*Bjt — -J- [{Mb* +Mb)Ai{0) + {Mb* — Mb)A2{0)\ . (14) 
Jit 

The Ai term is expected to dominate because the relative con¬ 
tribution of A 2 is suppressed by the ratio (Mg* — Mg)/(Mg* -f 
Mg) whose value is 0.004 for the physical B and B* masses. 
It is this relation that we use for our numerical calculation. 

















3 


TABLE I. Lattice simulation parameters. All ensembles are gen¬ 
erated using 2+1 flavours of domain-wall fermions and the Iwasaki 
gauge action. All valence pion masses are equal to the sea-pion mass. 


Lja 

aj fm 

mia 

rrisa 

# configs 

# sources 

M^/MeV 

24 

0.11 

0.005 

0.04 

1636 

1 

329 

24 

0.11 

0.010 

0.04 

1419 

1 

419 

24 

0.11 

0.020 

0.04 

345 

1 

558 

32 

0.08 

0.004 

0.03 

628 

2 

289 

32 

0.08 

0.006 

0.03 

889 

2 

345 

32 

0.08 

0.008 

0.03 

544 

2 

394 


III. CALCULATIONAL STRATEGY 

A. Light quarks and gauge fields 

Our analysis is carried out using ensembles produced by the 
RBC and UKQCD collaborations [25, 26] with the Iwasaki 
gauge action [27, 28] and 2-1-1 flavour dynamical domain- 
wall fermions [29, 30]. The configurations are at two lat¬ 
tice spacings, the finer 32^ ensembles have an inverse lattice 
spacing of = 2.281 (28) GeV and the coarser 24^ ensem¬ 
bles have = 1.729(25)GeV, corresponding to approxi¬ 

mately 0.086fm and 0.11 fm respectively. All ensembles have 
a spatial extent of 2.6fm. We simulate with unitary light- 
quarks corresponding to pion masses down to M-^ = 289MeV. 
On all ensembles the strange sea-quark mass is tuned to 
within 10% of its physical value. The fifth dimensional ex¬ 
tent of both lattices is Ls = 16, corresponding to a residual 
quark mass of = 0.003152 on the 24^ ensembles and 
mresfl = 0.0006664 on the 32^ ensembles. The lattice quark 
masses corresponding to the physical u/d and s quarks are 


=0.00136(4), ihsO = 0.0379(11) on the 24^ ensembles 
and iniiidQ = 0.00102(5), nisU = 0.0280(7) on the 32^ ensem¬ 
bles. The tildes indicate that these values include the residual 
quark mass. 

Our calculation makes use of unitary light-quark propaga¬ 
tors with point sources previously generated as part of the 
RBC/UKQCD Z?-physics program [31-35]. Full details of the 
ensembles and propagators used are presented in Table I. We 
perform a random translation on each gauge held configura¬ 
tion to minimise the effects of autocorrelations on our results, 
allowing us to use more closely spaced trajectories and gain 
statistics. For each configuration in the 32^ ensembles we use 
propagators computed at two time sources separated by half 
the lattice temporal extent to compensate for the smaller en¬ 
semble sizes. 


B. Bottom quarks 

Simulating heavy quarks on the lattice presents the problem 
of dealing with mga > 1. A number of approaches have been 
developed to tackle this problem. In the limit of infinite mass 
the quarks become a static source of colour charge and their 
lattice propagator reduces to a trace of a product of temporal 
gauge links. This is the static action of Eichten and Hill [36] 
which has been used extensively to calculate the coupling g„, 
most recently in [13, 37]. Another approach is nonrelativistic 
QCD [38], where the usual QCD Lagrangian is expanded in 
powers of v/c. 

Here we use the relativistic heavy-quark (RHQ) action [39- 
41] to simulate fully relativistic bottom quarks while con¬ 
trolling discretization effects. Although mga is large for the 
heavy quark in a heavy-light meson, the spatial momentum 
\pa\ is of G(aAQCD)- The RHQ action is an anisotropic Wil¬ 
son action with a Sheikholeslami-Wohlert term [42] 


SRHQ=a‘^Y.'^iy)[mo + YoDo + C7■D-^{Do)^-^CiD)^+Y.J^P<y^lvF^v) W(x). (15) 

x,y \ jiv J yyc 


El-KJiadra, Kronfeld and Mackenzie showed that, for cor¬ 
rectly tuned parameters, the anisotropic Clover action can be 
used to describe heavy quarks with controlled cut-off effects 
to all orders in mga and of 0{\pa\) [39]. Christ, Li, and 
Lin [41] later showed that only three independent parameters 
need to be determined and, further, presented a method for 
performing this parameter tuning non-perturbatively [43]. 

This tuning has been completed for b quarks [31] on the 
RBC/UKQCD configurations and those results (Table II) are 
exploited in this calculation. The heavy-quark propagators 
and the correlation functions used in this analysis are cal¬ 
culated using the CHROMA software library [44]. We apply 
Gaussian smearing to the heavy-quark propagators using pa¬ 
rameters tuned in [31]. Because the correlators become very 
small at large time separations owing to the large masses in the 


exponentials, we run the inverter for the heavy quark propa¬ 
gators to a very small relative residual (10^^^). We found that 
pursuing the conjugate gradient iteration to this small resid¬ 
ual is equivalent to demanding convergence separately for the 
residual on each time slice. 


C. Three-point correlation functions 

The matrix element that we wish to calculate in Eq. (12) 
corresponds to the quark-flow diagram shown in Eig. 1. To 
fully benefit from the available pre-calculated light-quark 
propagators, we arrange the calculation so that the axial- 
vector current is positioned at the light-quark propagator’s 
source. This means that we use the periodicity of the lattice. 
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TABLE II. Tuned RHQ parameters for b quarks from Ref. [31]. 
The uncertainties shown are statistical, heavy-quark discretization 
effects, lattice-scale uncertainty, and from experimental inputs (the 
spin-averaged -meson mass and 6, hyperfine splitting) respec¬ 
tively. 


fl/fm 

rtiQa 

Cp 

C 

0.11 

8.45(6)(13)(50)(7) 

5.8(1)(4)(4)(2) 

3.10(7)(11)(9)(0) 

0.08 

3.99(3)(6)(18)(3) 3.57(7)(22)(19)(14) 

1.93(4)(7)(3)(0) 


z = (0,0) 



FIG. 1. Quark flow diagram 

creating the B* meson at large time tx, propagating to the ori¬ 
gin at f = 0 , where the current acts and continuing to small 
time ty, where the Z?-meson is annihilated. With this setup, 
we need only one inversion, a heavy-quark propagator calcu¬ 
lated from ty, using a light-quark propagator for the sequential 
source. The necessary trace is 

L^^''^'"''^'"Tr[5,(0,y)755,(y,x)7v5,(^,0)7^75] . (16) 

x,y 

Because of the periodicity of the lattice there are further pos¬ 
sible contributions beyond the desired Wick contraction in 


Fig. 1 . An operator located at time t can also be considered at 
location t — T, where T is the lattice’s temporal extent. If we 
take into account separately the cases where > ty and < ty, 
there are eight possible contributing arrangements (consider¬ 
ing t + nT for integer n gives further contributions, but for 
increasing \n\ these are more and more suppressed). These 
are shown in Fig. 2. Using the approximation that the ground 
states immediately dominate the time dependence and that the 
matrix elements are all unity, we can estimate the time depen¬ 
dence of the three-point correlation function: 

d^\tx,ty) = 

Tty) “f U(UTr) “f ^x^ty, ( 17 ) 

^ (U; ty) “f ^ (U; ty) “f 6r(u, ty) H (t^, ty) tx ty. 

It is expected that the signal from which to extract the B*Bn 
coupling will be seen at large tx (approaching T from below), 
coming from the contribution shown as C{tx,ty) in Fig. 2. 
Figure 3 shows the relative size of the different contribu¬ 
tions as a function of tx with all matrix elements set to unity. 
As anticipated, C is the dominant contribution in a region 
T/2 + ty <tx <T. This result appears steady for a range of 
the masses Mb,Mb* ,M.^. Plotting the sum of the contributions 
as a function of U (setting all matrix elements to unity) we see 
a peak at tx = ty and an overall cosh-like form shifted by ty as 
shown on the right in Fig. 4. On the left of Fig. 4 we show 
good agreement with this form for our numerical data, with 
ty the time-position of the source for the sequential inversion. 
This gives us conhdence that we can extract the desired matrix 
element from the large U region of the three-point correlator. 


D. Ratios 

To access the matrix element in Eq. (12) we calculate the 
lattice three-point function 


^fiv {txity',p,p ) — 


T-tx>0 

-‘P'-HBiy)Ay{0)B;{x)) 


x,y 


1 / 21/2 

2Eb2Eb* 


{B{p')\Ay\B*ip,X))el^>e 


W„-EBty-Es*{T-tx] 


(18) 


and the vector and pseudoscalar two-point functions 


C^B^ {f,p)='£e-‘>’-^{B{x)B{0))^ZB — , 


Zb. 


g-Eg*t 

2Eb* 



PliPv \ 

7 ■ 


(19) 


( 20 ) 


If we set both the vector and pseudoscalar momenta to zero 
in Eq. (18), such that q = p = p' = 0 and = q^ = [Mb* — 


Mb)^ ~ 0, we can see from Eq. (12) that the only form factor 
accessible is Ai. Hence we form the ratio (not summed on i): 


f^l {tx,ty) 


cf{tx,ty-p = Q,p'=Q)z]l^z]l,^ 
4B(fv;F = 0)4?^,(r-u;p = 0) 
[Mb* +MB)Ai(q'o), 


( 21 ) 


where we can average over the three spatial directions (/ = 
1,2,3). To access the other form factors we need to inject 
a unit of momentum, such that q = p = ( 1 , 0 , 0 ) x 2%lL and 
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(X exp(-Mj;t^}exp ( - Mg.(T-tj) 


D{t^, fj,) oc exp ( - mJT-tj) exp ( - Mg(t^-t^)) 




G(t^, rp oc exp ( - Mg(T-t^)) exp(-Mg.r^) H{t^, t^) oc exp ( - m^(T-t^)) exp ( - 


FICj. 2. Lc^din^ contributions to the thi'ee“point correlutoi" for ty und <ty{E,F,G,H). 


p — 0. Following [6], we define the ratios: 


These allow access to the form factor A 2 through 




^3 


^4 


d^^iG,ty-p^Q,p' = Q)zf4i^ 

Y,{B{p')\Ao\B*ip,X))e[^\ (22) 

X 

cf^{t,,ty-p^Q,p' = Q)z]l^zf 
C^Bkh-^p’=0)d^s.{T-ta^Q) 
Y^{B[p')\Ai\B*{p,X))e^^\ (23) 

A 

dk{G,ty-p^Q,p' = Q)z]l^zf 

ck{ty-p' = o)c^}^.{T-t,-,pyo) 

j:{Bip')\A2\B*ip^))ek^ 

X 

{Mb* + Mb )Ai{q^). (24) 


A2 

Ai 


(Mb* +Mb)^ 

2Mlq] 


— q\+ Eb* {Eb* — Mb ) 


mI.{Eb*-Mb) R 2 Ml^qi Ri 
Eb* Ra Eb* Ra 


(25) 


The ratio in Eq. (25) is obtained at non-zero values of q^ and 
needs to be extrapolated to q^ = 0. However, from Eq. (14) 
the contribution of A 2 (0) relative to A 1 (0) is suppressed by the 
ratio {Mb* —Mb) / {Mb* XMb). The form factor Ai is obtained 
at q^ = {Mb* — Mb)^ from Eq. (21), but examination shows 
that the slight extrapolation to = 0 is not necessary at the 
resolution possible with the available statistics. If we define 
functions G\ and G 2 


Gi{q^) = {MB*+MB)Ai{q^), 

G2{q^) = {Mb* —MB)A2{q^), 
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FIG. 3. The estimated relative sizes, with ty = 8, of the eight contri- 
hutions to the B*Bn three-point function arising from different Wick 
contractions, using a linear scale (top) and log scale (bottom). Each 
contribution is scaled by dividing by the maximum contribution at 
that time. The matrix element of interest comes from C which is 
shown as a solid red line on both plots. This contribution dominates 
for large 


we can write the coupling as G\ (0) plus a small correction 
from the ratio G 2 /G 1 , 


gb 


Za 

2Mb 


Gi(0) 



(27) 


where Za is the light axial-vector current renormalization fac¬ 
tor. In our simulations A 2 is of comparable size to A]. The 
mass suppression in the ratio G 2 /G 1 means that the correc¬ 
tion term in (27) is at most 2% on our ensembles and typically 
at the sub percent level, with an error comparable to its size. 

We take Za from the RBC/UKQCD combined analysis of 
the light hadron spectrum, pseudoscalar meson decay con¬ 
stants and quark masses on the 24^ and 32^ ensembles [26]. 
The values are calculated from the ratio of the conserved and 
local vector currents, extrapolated to the chiral limit and are 
shown in Table III. 


TABLE III. Axial current renormalization factors used in this work, 
calculated in [26]. 


Ensemble 

a/fm 

Za 

24^ 

0.11 

0.7019(26) 

32^ 

0.086 

0.7396(17) 


TABLE IV. Eit ranges used for the two-point functions and the ra¬ 
tios. Eor non-zero momenta, equivalent combinations are averaged. 
Eit ranges are the same for different light-quark masses at the same 
lattice spacing, except in the case of the B* two-point function with 
momentum (1,1,0) on 24^, m/a = 0.020 for which a range 7-15 is 
used in place of 9-15 listed in the table. 



pa/{2n/L) 

fit range 

243 

^max 

— 223 - 

B 

(0,0,0) 

8-16 

8-17 

B* 

(0,0,0) 

7-15 

8-16 

Rl 

(0,0,0) 

50-58 

50-58 

B* 

(1,0,0) 

7-15 

8-16 

Rl 

(1,0,0) 

51-60 

47-55 

Rl 

(1,0,0) 

50-56 

46-55 

R 4 

(1,0,0) 

50-57 

47-56 

B* 

(1,1,0) 

9-15 

10-16 

Rl 

(1,1,0) 

51-60 

47-55 

Rl 

(1,1,0) 

49-55 

46-55 

R 4 

(1,1,0) 

51-57 

46-54 


IV. ANALYSIS 


A. Correlator fits 

We first calculate the three-point function on the 24^ en¬ 
semble with anil = 0.005 for values of ty ranging from 6 to 
18. Examining the data for R\ (see Fig. 5), it is clear that the 
best signal is achieved with ty = 6. We therefore choose ty = 6 
for our analysis on the 24^ ensembles and ty = 8 on the 32^ en¬ 
sembles because it corresponds to the same physical distance. 

Fig. 6 shows the ratios Ri, R 2 , R 3 and R 4 calculated on the 
243 , m/a = 0.005 ensemble, and Fig. 7 shows the vector and 
pseudoscalar effective-mass plots. In all cases, we estimate 
the statistical error with a single-elimination jackknife. The 
two-point functions are fitted to a single exponential to extract 
Zb and Zb*. Using these values of Z^(*), we then fit the ratios 
to a constant in the regions given in Table IV where we expect 
excited-state contamination to be small. We choose the fit 
ranges for each ratio such that we obtain a good correlated 
X^/dof, and apply them to all ensembles of the same lattice 
spacing consistently. 

From the ratios we use the procedure described in the previ¬ 
ous section to extract gb on each ensemble, giving the values 
listed in Table V. 
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FIG. 4. The three-point correlator of Eq. (18) evaluated on the 24^, 
matches that predicted from the analysis in Sec. Ill C (right). 



FIG. 5. The ratio Ri {t,ty) evaluated for different values of ty on the 
24^, m/fl = 0.005 ensemble, ty = 6 gives the cleanest signal and 
longest plateau. The points for different ty have small horizontal off¬ 
sets to help distinguish them on the plot. 


TABLE V. Results for gi, on the 24^ and 32^ ensembles. Errors for 
gl, are statistical. Pion masses are from [26]. 


Lja 

mia 

M^/GeV^ 

gh 

24 

0.005 

0.108 

0.533 ±0.027 

24 

0.010 

0.175 

0.568 ±0.023 

24 

0.020 

0.311 

0.580 ±0.035 

32 

0.004 

0.084 

0.548 ±0.020 

32 

0.006 

0.119 

0.603 ±0.016 

32 

0.008 

0.155 

0.596±0.018 


= 0.005 ensemble with ty = 6 (left). The time dependence closely 


B. Chiral and Continuum Extrapolations 

We perform a chiral extrapolation using the SU(2) HM;kPT 
formula for the axial coupling matrix element derived in [45]; 

+ (28) 

which is next-to-leading order in the chiral expansion, but 
only leading order in the heavy-quark expansion. We param¬ 
eterize the light-quark and gluon discretization effects with 
an term, as expected for the domain-wall light-quark and 
Iwasaki gauge actions. The lattice-spacing dependence from 
the RHQ action is more complicated. However, we argue 
in the next section that heavy-quark discretization effects are 
negligible and that an extrapolation in captures the lead¬ 
ing scaling behaviour. We use the PDG value [46] of the pion 
decay constant, fn = 130.4MeV. 

Fig. 8 shows the chiral-continuum extrapolation of the nu¬ 
merical simulation data to the physical light-quark mass and 
continuum using Eq. (28). The values of Za, as calculated 
in [25] and [26], are included for each ensemble. The statis¬ 
tical errors in Za are added in quadrature to the Monte-Carlo 
statistical errors in the lattice data for gh before performing 
the chiral fit. The six ensembles are statistically independent, 
hence the fit is calculated by minimizing an uncorrelated 
function. The parameters that minimise the x^ are go = 0.515, 
a = -1.324GeV^^ j3 = -0.648GeVl We use these fitted 
parameters and their covariance matrix to estimate the error 
bands in our plots and to give the errors in gy. Our fitted result 
is gh = 0.557 ± 0.027. 

V. SYSTEMATIC ERRORS 
A. Chiral Extrapolation 


Our chiral extrapolation relies on NLO SU(2) HMjPT with 
pion masses above 400MeV. Therefore we may be using the 
theory beyond its range of applicability and we are certainly 
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FIG. 6. Ratios R\{t,ty), R 2 {t,ty), R^{t,ty) and Ri\[t,ty) for fy = 6 on the 24^, mia = 0.005 ensemble. Ri is calculated with p = p' = 0, the other 
ratios with p = (l,0,0,0)2;r/L, p' = 0. 




FIG. 7. 6-meson (left) and 6*-meson (right) effective masses on the 24^, m/fl = 0.005 ensemble. 


omitting higher order terms in the chiral expansion. To es¬ 
timate the uncertainty this introduces we consider a range of 
possible fits. First, we consider the effect of neglecting the 
heaviest mass from each ensemble (center left plot in Fig. 9). 
This alters the form of the fit dramatically but does not sig¬ 
nificantly change the final result. In the bottom row of Fig. 9 
we replace in the coefficient of the NLO chiral logarithms 
with /^ = 156.1 MeV [46] or with /o = 112MeV in the SU(2) 


chiral limit from the RBC/UKQCD light pseudoscalar me¬ 
son analysis [26]. This changes the relative size of NLO and 
NNLO and higher-order terms in the chiral expansion. Fi¬ 
nally, we note that our data does not show any strong evidence 
of chiral log curvature, presumably because our lightest data 
point corresponds to Mjc « 289 MeV and is still rather heavy. 
We therefore consider an analytic fit, shown in the centre right 
plot of Fig. 9, where we extrapolate linearly in M^. Of these 
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ment in partially quenched HM;ifPT. This allows us to eval¬ 
uate the expression with different valence and sea strange- 
quark masses. The NLO matrix element has four different 
contributions, coming from so called sunset diagrams, wave- 
function renormalization, tadpole diagrams and the NLO an¬ 
alytic terms. We have calculated the effect of a 10% change 
in the sea strange-quark mass in the loop diagrams, assum¬ 
ing the values of the low-energy constants obtained from our 
preferred chiral fit, on the value of the coupling gi,. We find 
a change in gt, of 1.5%. This result is numerically consis¬ 
tent with the effect of the strange sea-quark mass on the pion 
decay constant observed by the RBC/UKQCD collaboration 
in [26]. Therefore we ascribe an error of 1.5% in gh due to the 
unphysical strange-quark mass. 


FIG. 8. Chiral and continuum extrapolation. The bottom (blue) curve 
is the fit through the 24^ ensemble points. The (red) curve above is 
the fit through the 32^ ensemble points and the (green) error band and 
curve show the continuum extrapolation. The intersection with the 
vertical grey line corresponds to the physical pion mass. The result 
at the physical mass is shown by the (black) point with error. 


variations, the largest difference from our central value for 
gh is from the linear fit in and a^. This value is larger 
than our full chiral-continuum fit by 10.6%. Because the chi¬ 
ral and continuum extrapolations are treated together in our 
fitting procedure, however, discretization and chiral extrapo¬ 
lation errors cannot be fully disentangled. In section V E 2 be¬ 
low we consider light-quark and gluon discretization errors, 
estimating a systematic error of 11.5%. This is the largest de¬ 
viation seen in the chiral-continuum fits in Fig. 9 and is there¬ 
fore the error we take for the combined chiral and continuum 
extrapolation. 


B. Lattice-scale uncertainty 

The coupling gh is a dimensionless number calculated from 
ratios of correlators, so it should have only a mild dependence 
on the physical value of the lattice spacing. However, varia¬ 
tions in a affect the chiral and continuum extrapolations. We 
estimate the error in gh due to the lattice-spacing uncertainty 
by varying the 24^ and 32^ lattice spacings by their quoted 
(statistical plus systematic) uncertainties, O 24 and 032 [26], 
one at a time whilst keeping the other fixed. Shifting the lat¬ 
tice spacing on the finest ensemble changes gh by 0.7%, and 
on the coarser ensemble gh changes by 0.6%. Therefore as¬ 
cribing an error of 0.9% (the sum in quadrature) to this source 
of uncertainty seems a conservative estimate. 


D. RHQ parameter uncertainties 

1. Statistical 

To test the dependence of gh on the uncertainties in the 
tuned RHQ parameters we calculate the coupling on the 24^ 
m/a — 0.005 ensemble using the full “box” of RHQ parame¬ 
ters used to interpolate to the tuned values: 


moa 


ttlQQ zb UfyiQa 


moa 


moa 

Cp 

7 

Cp 

7 

Cp ± (Tcp 

7 

Cp 

c . 


c 


c 




(29) 

For our 24^ ensembles, the box parameters are given by 

{moa,Cp,Q = (8.40,5.80,3.20), 

= (0.15,0.45,0.30). 

We then linearly interpolate gh to the point of the tuned pa¬ 
rameters. By following this procedure underneath the jack¬ 
knife we can propagate the statistical errors from parameter 
tuning through to gh- Comparison of this determination to the 
result calculated directly using the tuned values of the parame¬ 
ters gives a measure of how sensitive gh is to the uncertainties 
arising from the tuning. We find that the central values differ 
by 0.01% and the errors agree to two significant figures. In 
the context of the overall uncertainty this can be considered 
negligible. 

Figure 10 shows gh calculated on the seven sets of param¬ 
eters indicated in Eq. (29) for the 24^ m/a = 0.005 ensemble. 


2. Systematic 


C. Unphysical sea strange-quark mass 

Our simulation is performed with a sea strange-quark mass 
that differs from the physical value by approximately 10%. To 
investigate the effect of the sea strange-quark mass on gh we 
use results from [45] for the NLO axial current matrix ele- 


We also consider the effect on gh of systematic uncertainties 
in the RHQ parameters. These are estimated in Ref. [31] and 
given in Table IT The three significant contributors are heavy- 
quark discretization effects, uncertainty in the lattice spacing, 
and uncertainty from the experimental inputs. To determine 
the sensitivity of gh to these uncertainties we use the calcu¬ 
lation on the box of parameters, Eq. (29), described in the 
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mI /GeV^ mI /GeV^ 


FIG. 9. Variations of chiral fits. Top left: 32^ data points only; top right: data points from both lattices with no a? term. Center left: heaviest 
masses dropped; center right: fit to a function linear in M^. Bottom left: replacing fjc with = 156.1MeV; bottom right: replacing fn with 
/o = 112MeV. In each plot, the result of the preferred fit from Fig. 8 is shown as the black point, with error, on the vertical line at the physical 
pion mass. 


previous subsection. We assume a linear dependence of gi, on 
the RHQ parameters for small shifts, then shift one parameter 
at a time by each systematic uncertainty and take the overall 
error as the effect of each of these shifts added in quadrature. 
The combined effect, shown in Table VI, is an error of 1.5% 
in gb- 


E. Discretization errors 

1. Heavy-quark discretization errors 

We estimate heavy-quark discretization errors using an ef¬ 
fective field theory approach [39, 47, 48] in which both our 



















11 


0.57 

0.56 

0.55 

0.54 

'o.53 

0.52 

0.51 

0.50 

0.57 

0.56 

0.55 

0.54 

'o.53 

0.52 

0.51 

0.50 

0.57 

0.56 

0.55 

0.54 

'o.53 

0.52 

0.51 

0.50 


8.25 8.30 8.35 8.40 8.45 8.50 8.55 

mo a 


5.4 


5.6 


5.8 


6.0 


6.2 


2.9 3.0 3.1 


3.2 


3.3 3.4 


3.5 


FIG. 10. gh calculated for the sets of RHQ parameters used to define 
the parameter “box” on the 24^ m;a = 0.005 ensemble. The blue 
points are the results for the box parameter choices and the red point 
shows gj, calculated directly at the tuned parameter values. 


lattice theory and QCD are described by effective continuum 
Lagrangians built from the same operators and errors stem 
from mismatches between the short-distance coefficients of 
higher-dimension operators in the two effective theories. Ok¬ 
tay and Kronfeld [48] have catalogued the relevant operators 
and calculated the mismatch coefficients at tree level. 

Because we are evaluating a matrix element of the light- 


TABLE VI. The effect of systematic uncertainties in the RHQ pa¬ 
rameters on gh- Each parameter was shifted by the uncertainty from 
each source and the effect on gi, calculated by assuming gi, depends 
linearly on the parameters. 


Source 


% error from parameter 

moa 

Cp 

c 

Total 

HQ discretization 

0.25 

0.65 

0.30 

0.76 

Lattice scale 

0.97 

0.65 

0.24 

1.19 

Experimental inputs 

0.14 

0.33 

0 

0.35 

Total 

1.01 

0.98 

0.38 

1.46 


quark axial-vector current, heavy-quark discretization errors 
stem from mismatches in higher-dimension operators in the 
heavy-quark action which correct the B and B* meson masses. 
We expect these effects to be negligible. From our tuning 
procedure [31] we can relate changes in the meson masses 
to changes in the RHQ parameters moa, Cp and i^, while in 
section V D below, we relate changes in the RHQ parameters 
to changes in gy. Hence we can estimate the effect of errors in 
the meson masses on gi,. 

In Appendix C of [31], we estimated the heavy-quark 
discretization error on the spin-averaged Bj meson mass as 
0.05%. Also in [31], that spin-averaged mass was most sen¬ 
sitive to variations in mpa, with a 0.05% shift corresponding 
to a change of around 0.02 in moa on the 24^ m/a = 0.005 
ensemble. From section V D, shifting moa by the halfwidth 
of our tuning “box” changes gh by no more than 1.5%. For 
the 24^ m/a = 0.005 ensemble, this shift in moa is 0.15 and 
hence we expect a heavy-quark discretization error on gy of 
no more than (0.02/0.15) x 1.5% = 0.2%, which is negligi¬ 
ble compared to our overall uncertainty. 


2. Light-quark and gluon discretization errors 

Leading discretization errors from the domain-wall light- 
quark action and the Iwasaki gauge action are both 0{a^) and 
are included as an term in the combined chiral-continuum 
extrapolation. However the data is also compatible within er¬ 
rors with assuming no lattice-spacing dependence; a fit with 
no term also yields an acceptable, albeit larger, j^/dof. 
The top row of Fig. 9 shows chiral fits to the data without 
an term. To estimate the systematic errors coming from the 
continuum extrapolation we use the difference of 11.5% in gt, 
between a fit to our finest data set (a « 0.086fm) and the 
extrapolation using both lattice spacings. This is the largest 
effect in all variations of our chiral and continuum extrapola¬ 
tions and is therefore the value appearing in Table VII for the 
combined chiral and continuum extrapolation uncertainty. 


F. Finite-volume effects 

We expect that finite-volume effects are small since there 
are no propagating light particles in the simulated system. To 
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TABLE VII. Error budget for systematic and statistical errors. 


Statistical errors 

4.8% 

Chiral and continuum extrapolation 

11.5% 

Lattice scale uncertainty 

0.9% 

Finite volume effects 

1 .0% 

RHQ parameter uncertainties 

1.5% 

Unphysical sea strange-quark mass 

1.5% 

Systematic errors total 

11 .8% 


charm — 

Becirevic and Sanfilippo 2012 [8] 

— Becirevic and Haas 2011 [49] 

— Abada et al 2002 [6] 

RBC-UKQCD this work 

bottom 

Abada et al 2004 (interpolation) [7] 

, 

Detmold et al 2012 [12] 


ALPHA 2014 [13] 

* 

ALPHA 2010 [50] 

Becirevic et al 2009 [51] 

static —*— 

Ohki etal 2008 [11] 


Negishi et al 2007 [52] 

Abada et al 2004 [7] 

UKQCD 1998 [10] 


03 04 O'S 0^6 Q.l O'S 


FIG. 11. Comparison of quenched (downward triangles), Nf = 2 (up¬ 
ward triangles) and Nf = 2+1 (circles) calculations of the couplings 
gc, gh and g [6-8, 10-13, 49-52]. Error bars represent the sum in 
quadrature of all quoted errors (statistical and systematic). 

estimate their size we compare the value of gi, obtained from 
an NLO heavy-meson fit to our data, with and without 
hnite volume effects included. We compare the finite and 
inhnite-volume ht result at all of our simulated pion-mass val¬ 
ues. The largest finite-volume correction, which occurs for 
our lightest pion mass, is ^ 1%, so we take 1% as the finite- 
volume error in our calculation of gt,. 


VI. CONCLUSIONS 

The sum in quadrature of all the systematic errors described 
in section V gives a total systematic uncertainty of 12%. Our 
final error budget is given in Table VII and our final value of 
the coupling gi, including statistical and systematic errors is 

8 b = 0.56(3) stat (7) sys (31) 

Our calculation is the first directly at the physical h-quark 
mass, and has a complete systematic error budget. Fig. 11 
compares our result to earlier dynamical calculations at the 
charm-quark mass and in the static limit. The dependence 
of g on the value of the heavy-quark mass is mild, and our 
result lies in the region that would be expected from interpo¬ 
lating between the charm- and infinite-mass determinations. 
Our result is compatible with the experimental value = 
0.570 ± 0.004 ± 0.005 extracted from the natural linewidth 
of the transition D*(2010)+ -> by the BaBar Collab¬ 
oration in [17]. This further suggests that l/wg corrections 
to the coupling g are small. Our result has been used by 
the RBC/UKQCD collaboration in the chiral extrapolations 
of numerical lattice data for the B-meson leptonic decay con¬ 
stants [32, 33] and B -+ niv and Bs -+ Kiv semileptonic form 
factors [34, 35]. 
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